Magnetization process of spin ice in a [111] magnetic field 
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Spin ice in a magnetic field in the [111] direction displays two magnetization plateaux, one at 
saturation and an intermediate one with finite entropy. We study the crossovers between the different 
regimes from a point of view of (entropically) interacting defects. We develop an analytical theory 
for the nearest-neighbor spin ice model, which covers most of the magnetization curve. We find 
that the entropy is non-monotonic, exhibiting a giant spike between the two plateaux. This regime 
is described by a monomer-dimer model with tunable fugacities. At low fields, we develop an RG 
treatment for the extended string defects, and we compare our results to extensive Monte Carlo 
simulations. We address the implications of our results for cooling by adiabatic (de)magnetization. 
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I. INTRODUCTION 

Recent experiments on the spin ice compounds^ 
Ho2Ti2C>7 and Dy2Ti2C<7 have uncovered an intriguing 
set of phenomena when unicrystalline samples are placed 
in an external magnetic field in the [111] directioniiii^& 
For a review on spin ice, see Ref. |jj 

The discovery of a plateau in the magnetization be- 
low saturation, first predicted theoretically and explored 
in Monte Carlo simulations^^ has been particularly re- 
markable as it was found to retain a fraction of the zero- 
field spin ice entropyAifiaii In this regime, the system 
is well described by a two-dimensional Ising model on a 
kagome lattice in a longitudinal field, which is in turn 
equivalent to a hexagonal lattice dimer model»i2*ii*i£. 

Recently, two of the present authors have studied the 
thermodynamics and correlations of the [111] plateau. 11 
This work has led to the identification of mechanisms 
which terminate the plateau. At the high-field end, the 
termination occurs via the proliferation of monomer de- 
fects in the underlying dimer model. At low fields, a more 
exotic extended string defect restores three dimensional- 
ity. The asymptotic density of both kinds of defects was 
estimated in Ref. 

In this paper, we consider in detail the full magnetiza- 
tion curve from zero-field to saturation. A brief synopsis 
of the exotic thermodynamic properties of spin ice is in 
a [111] field is sketched in Fig.^ The aim of this paper 
is to identify the different regimes of the magnetization 
curves, to provide analytical theories for them, and to 
test them against Monte Carlo simulations, and finally 
against experiment. 

Near zero field, we use the accurate self-consistent 
Hartree approximation^ to provide an analytical approx- 
imation for the linear response regime. At the low field 
end of the plateau, we develop mean field and renormal- 
ization group treatments for the extended string defects, 
which we use to analyze the in-plane and out-of-plane 
correlations. We compare these with Monte Carlo sim- 
ulations using an efficient cluster algorithm, which al- 
lows us to obtain accurate data from the linear response 
regime to the beginning of the [111] plateau. We find 



that the mean field treatment is accurate at the lowest 
fields, where the string density would be relatively high. 
The renormalization group treatment compares well with 
simulation in the dilute string limit. At even higher fields, 
the plateau is approached and the suppression of the en- 
tropic activation of strings becomes apparent as a finite- 
size effect. 

At the high-field termination of the plateau, we ob- 
serve a giant peak in the entropy, which even exceeds the 
zero field Pauling value, despite the fact that a quarter 
of all spins are pinned. We model this phenomenon by 
a monomer-dimer model on the honeycomb lattice with 
varying fugacities. (At the point where the all fugacities 
equal 1, this model turns out to be one of 'hard bow-ties' 
on the kagome lattice.) We analyze this model within 
a Bcthe approximation and also by using results from a 
high-order series expansion^ 

We show that the entropy peak is due to the crossing of 
an extensive number of energy levels which have macro- 
scopic entropies. Comparing this theory with Monte 
Carlo simulations of the appropriate monomer-dimer 
model, we find that the simple Bethe approximation is 
accurate for moderate to large monomer densities. 

We point out that this theory predicts to a crossing 
point in the plots of magnetisation versus field at differ- 
ent temperatures. In addition, there is a further cross- 
ing point at lower fields, where the corrections to the 
magnetisation due to monomer and string defects almost 
cancel one another. 

We then address the connection of these results to ex- 
periment, in particular pointing out the presence of (at 
least a vestige) of the entropy peak in existing experi- 
mental data. 

We then discuss the implications of the entropy peak 
for magnetocaloric manipulations. In particular, we ar- 
gue that it arises in a more general set of models. It can, 
in principle, be used to effect cooling in a field, both by 
adiabatic demagnetization, and by adiabatic magnetiza- 
tion. Finally, we close with some concluding remarks. 
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FIG. 1: Properties of spin-ice as the [111] magnetic field is 
varied. These curves are for illustration and do not show 
actual numerical or experimental data. We have indicated 
the regions where various analytic approaches discussed in 
the text apply. 

II. MODEL AND NOTATION 

A general model of spin ice includes the single-ion 
anisotropy, the exchange interaction, and the dipolar in- 
teraction. In this work we use a simplified modei in 
which the long-range dipolar interaction is truncated be- 
yond the nearest-neighbor spins. While the exchange 
interaction in spin ice compounds is antiferromagnetic, 
the effective interaction (exchange plus nearest-neighbor 
dipolar) is ferromagnetic. The Hamiltonian for unit- 
length spins Sj may be written as 

<ij> i 

+ g/iflJ^B-Si, (2.1) 

i 

where J' oS is an effective exchange coupling. The sec- 
ond term is the easy axis anisotropy of strength E < 0, 




FIG. 2: The pyrochlore lattice of corner-sharing tetrahedra. 




FIG. 3: A single tetrahedron inscribed in a cube. The easy 
axes of the pyrochlore lattice (or (111) axes), d K , are indicated 
by the short-dashed lines. 



\E\ > 50 A" , which is much larger than the exchange and 
dipolar interaction strengths. The unit vectors d K (i) are 
the local easy axes of the pyrochlore lattice (see Fig. |2 
and Fig. The third term is the interaction with a 
magnetic field of strength B, g/J-sJ being the magnetic 
moment of the spins. Both experiment and theory indi- 
cate that this simplified model is a good description of 
spin ice at moderate temperatures. 

In our analysis, we take the single ion anisotropy to be 
infinite so the spins are constrained to lie along their local 
easy axes. In this limit, it is convenient to describe the 
system by the Ising pseudospins di, where = <7jd re (i). 
The pseudospin a t =+l(-l) if the physical spin points 
into (out of) its associated up-pointing tetrahedron. We 
may write an effective Hamiltonian for the pseudospins: 

H = J cS ^2 (J i (T 3 ~ 9VbJ Y B • d K ( i) cr i , (2.2) 

<ij> i 

where J c g = J^ ff /3 > 0. 
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III. THE LOW FIELD REGIME 



thermodynamics 



At zero magnetic field and zero temperature, the ferro- 
magnetic interaction gives rise to an "ice rule" constraint: 
the pseudospins on each tetrahedron must sum to zero, 
Ere "*! = 0- I n t erms °f t ne physical spins, on each 
tetrahedron two will point inwards (towards the center) 
and two will point outwards (away from the center). The 
set of configurations satisfying the ice rule comprises the 
zero-field spin ice ground state manifold. At low mag- 
netic fields (and low temperatures), the system will con- 
tinue to obey the ice rule, though the magnetic field will 
favor certain states among those in the zero-field ground 
state manifold. 

We have performed extensive Monte Carlo simulations 
of the low field regime, from zero-field up till the low 
field plateau termination, using a loop algorithm, which 
is discussed in Appendix A. Our algorithm probes only 
spin ice ground states (two spins in and two out on each 
tetrahedron) and is thus applicable at low temperatures 
T <C J ff and low magnetic fields, where the density of 
monomer defects, which are responsible for the high field 
plateau termination, is low. The simulation is written 
in terms of a pyrochlore lattice with the conventional 16 
site cubic unit cell (which contains four tetrahedra of each 
kind). The simulations have been done for systems with 
16, 128, 432, 1024, 2000, 3456, 5488, 8192, and 16000 
sites. For a system with 16000 sites, we perform 2.5 x 10 6 
loop flips for equilibration and 5 x 10 7 for averaging. For 
other system sizes, we we perform 1 x 10 7 loop flips for 
equilibration and 2 x 10 8 for averaging. The simulated 
magnetization as a function of the magnetic field strength 
is shown in Fig.0] The magnetization attains the plateau 
value at fields much larger than the temperature. 



dU 

T ' T ' 
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(3.1) 
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FIG. 4: The magnetization from Monte Carlo simulations. 



A. The linear response regime 

We may calculate the ground state entropy of spin ice 
at zero field by numerically integrating the first law of 



Noting that the magnetization is constant and equal to 
— g/isJ/3 per spin on the plateau and is zero at zero 
field, and that the value of the entropy on the plateau 
is S/k B = 0.08076512*11, we obtain for the entropy of 
spin ice, S/ks — 0.2051 ± 0.0001. Our value is very close 
to Pauling's estimate S/ks = 0.202733 and is consis- 
tent with the most accurate current theoretical estimate 
S/k B = 0.20501 ± 0.00005-i^ 

At zero field, we use the self-consistent Hartree ap- 
proximation, which is known to give a quantitatively ac- 
curate approximation to the ground state correlations 
of spin ice»i£ This gives x — 2(g^s J) 2 /3/cbT for spin 
ice. This compares well with our Monte Carlo result, 
X = (0.66735 ± 0.0003)(c^ B J) 2 /£: B T for a system with 
16000 sites. 



B. String defects and their interactions 

1. General description 

Figure presents the underlying pyrochlore lattice of 
spin ice and fig. |21 shows the [111] direction. It is con- 
venient to visualize the pyrochlore lattice as a stack of 
alternating kagome and triangular planes, the [111] direc- 
tion being the direction in which the planes are stacked. 
Each spin lies on a corner shared by an up-pointing and 
down-pointing tetrahedron. 

If the [111] magnetic field is large enough, the spins 
in the triangular planes align with the field; the kagome 
planes decouple from one another; and the system is well 
described by a two-dimensional model. This describes 
spin ice on the plateau. At fields slightly lower than 
the plateau, excitations called string defects?*^ restore 
three-dimensionality and are responsible for the low field 
termination of the plateau. 

To describe these defects, we consider the entropic ben- 
efit of relaxing the condition that the triangular planes 
are polarized. Suppose we flip a spin in some triangu- 
lar layer. Then, by the ice rule constraint, we must also 
flip a spin in each of the two neighboring kagome layers 
(on the two tetrahedra that are sharing the first flipped 
spin) . Flipping these kagome spins requires flipping spins 
in each of the two neighboring triangular layers, which 
requires flipping spins in the two next-nearest kagome 
layers and so on. The resulting "string defect" is an ex- 
citation that extends through the system. The energy 
cost, per kagome-triangle bilayer, of creating the string is 
E s = Sg/iB JB/3. To estimate the entropy, we note that 
creating a string actually involves creating a pair of de- 
fects in each kagome plane. A "positive" defect connects 
the kagome plane to the kagome plane directly above it 
via a flipped spin in the intermediate triangular plane. 
Similarly, a "negative" defect connects the kagome plane 
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to the kagome plane directly below it. These two de- 
fects may be separated by flipping pairs of spins point- 
ing in different directions on neighboring triangles of the 
kagome plane. The entropy in the kagome plane depends 
on this separation, which is the basis for the interaction 
between defects discussed below. Ignoring this correc- 
tion, the positive defect may be placed anywhere in the 
plane (which fixes the position of the negative defect in 
the layer above). This implies that the entropy per bi- 
layer is S ~ In A, where A is the area of a layer. This 
shows that for a given magnetic field, string defects are 
favored in a sufficiently large system. For a given sys- 
tem size, strings are favored at sufficiently low magnetic 
fields. 



2. Interactions 



(spin-wave) and "monomer" (vortex) contributions. A 
standard calculation^ gives the entropy of the monomer 
piece: 



S m = ^ / / 



'a{r)a{?) C 
\J J d 2 rd 2 r'a(r)a{P)^-Kh^- 



r — r 

T 

f—r'\ 



(3.4) 

where k = 1/2 and r is a hard-core radius comparable 
to the lattice spacing. This shows that the entropic in- 
teraction between two defects separated by distance r is 
given by piP2V(\fl — r5|) where is +1 (-1) for a positive 
(negative) defect and V(R) = —k\u(R/t). 



For magnetic fields in the plateau region, the triangular 
spins are fixed while each kagome plane contains two up 
pseudospins (<r = 1) and one down pseudospin (a = —1). 
This Ising model on the kagome lattice may be mapped 
onto the dimer model on the hexagonal latticei2*ii*i£, 
where a down pseudospin corresponds to a dimer on the 
hexagonal lattice. In this language, a string defect ap- 
pears as a pair of oppos itely charged monomers. 

As discussed in Ref. Illl a monomer-dimer covering 
may be described by assigning a height variable hi to 
each site i of the triangular lattice dual to the hexagonal 
lattice on which the dimers lie. The heights are assigned 
as follows. Moving from a site to a nearest neighbor site 
by moving clockwise around an up- (down-) triangle will 
increase (decrease) the height by +2 (-2) if a dimer is 
crossed. If a dimer is not crossed, then the height will 
decrease (increase) by -1 (+1). According to these rules, 
traversing a closed loop in the dual lattice will result 
in a height difference of +3 (-3) if a positive (negative) 
monomer is enclosed and otherwise. We note that the 
overall sign of the height assignments is a matter of con- 
vention and we may as well have chosen the hi so that 
traversing a closed loop containing a positive (negative) 
monomer gave a height difference of -3. 

In a coarse-grained description, the hi are replace d by 
a real, continuum field h(f) and as discussed in Ref. [T^, 
the entropy associated with a height field h(r) is given to 
lowest order ingredients by: 



S 



where K = § for the honeycomb lattice 
field has the property: 
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(3.2) 



The height 



V/i • dr 



3 / d 2 ra(r) 
Is 



(3.3) 



where u{r) is the monomer charge density and S is the 
region enclosed by the loop C. We may proceed by anal- 
ogy with the 2d XY modelii and divide h into "dimer" 



3. Mean field calculation 

If the number of defects is fairly large, we may expect 
the interaction to be sufficiently screened to justify the 
use of variational mean field theory^ We will investigate 
the in-plane and out-of-plane correlations for the defects. 

We consider a layered system of two-dimensional 
planes (indexed by the label k which ranges from — K 
to K) where each plane contains N positive and N neg- 
ative defects (which we refer to as charges) that interact 
logarithmically. The string constraint requires that each 
positive charge in layer k is rigidly connected to a neg- 
ative charge in the layer k + 1. We impose a periodic 
boundary condition to connect the positive charges in 
the K th layer to the negative charges in the — ATth layer. 

We formally impose the constraint by writing the 
"Hamiltonian" in terms of positive charges alone. The 
planes are stacked in the z-direction. Let x\ be the in- 
plane position of the ith positive charge in the /cth layer. 
In absence of external fields, the entropy of a particular 
configuration of N defects is given by: 



K 



N 



H 



k=-K ijij 



N 

£ 



V(\x1-x) +1 \)) (3.5) 



Here V(R) — —k\yi(R/t), where r is a hard-core radius 
defining the minimum separation between two charges 
and k = 1/2. The first term corresponds to the repul- 
sion of positive charges within the same layer. The ab- 
sence of a factor of i in front of this term is due to the 
string constraint: bringing two positive charges in the 
same plane close together also involves bringing together 
their negative partners in the plane above. In terms of 
our positive charge formulation, this means the repulsion 
is twice as large. The second term is the interlayer in- 
teraction. Physically, a positive charge in layer k has a 
negative partner in the layer k + 1 which attracts the pos- 
itive charges in layer k+1. In terms of our positive charge 
formulation, charges repel charges in the same plane but 
attract charges in nearest neighbor planes. 
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We assume a variational mean field density of the form: 

P k (Xi) 



p{x 1 



K 



k=-Ki=\ 



n n 



N 



(3.6) 



which asserts that all particles in a given layer fc have the 
same probability density p k (x)/N, but the density may 
vary from layer to layer. We also need the normalizing 
condition: 



d z xp k {x) = N 



(3.7) 



This trial function implies a variational entropy func- 
tional: 

S p ,n = £ (-~ f f d 2 xd 2 x'(p k (x)- P k+1 (x)) 



k=-K 

x (p k (x')- p k+1 (x'))V(\x-x'\ 
dV(*)k(^)) 



(3.8) 



This functional is maximized when the density is uniform 
p k (x) = % which gives S PtN = (2K+l)N\nA. To inves- 
tigate the linear response of the system, we may apply a 
perturbing potential to the objects in the k — plane. In 
particular, we consider the effect on the density of placing 
a positive charge at the origin of the plane. The details 
are given in Ref. but we may quote the result: 



s 2 (s 2 + 2) 



l + s 2 + 2 N /s 2 (s 2 + 2) l + s 2 + v /s 2 (s 2 +2) 



fc-1 



(3.9) 



where the in-plane length scale is given by £|| = 

,1/2 

. We note first that this expression diverges 



V AttkN 

at small x for k = 0, which is not surprising because the 
assumption of a linear response would be not be valid so 
close to the perturbing charge. The expression would be 
valid at larger k and an interesting feature is that when 
x = the decay in the z-direction does not depend on 
any physical parameters, i.e. there is no length scale in 
the z direction. We will return to this point in the next 
section. 

To connect with our physical problem, we note that 
at a given temperature, we will have an expected value 
of defects which may be calculated from the partition 
function: 



Z 



-0A 



(2K+l)N 

E — e A 



jV 



(AH) 



\\2K+1 



(3.10) 



where Sn is the entropy of having N defects and y = 
e -E s /k B T j g f U g ac ity f a positive defect (y 2K+1 is the 



JV 

as 



fugacity of a "string"). In mean field, we may replace S 
by S p ,n = (2K + 1)./Vln A From this, we may shoi 
that < N >^ yA and using our earlier expression, we 
find that: 



\\,MF 



e?q>(8gpB JB/3ksT) 



(3.11) 



4- RG calculation 



When the gas of defects is fairly dilute, we may ex- 
pect that the screening is not effective enough to jus- 
tify a mean field treatment. In this section, we account 
for fluctuations by making a real space renormalization 
group calculation using methods similar to the Kosterlitz 
treatment of the 2d Coulomb gas»i2^ 

The dynamical objects described by Hamiltonian 13.51 
are dipoles of length 1 . We need to generalize this model 
in order to do an RG calculation. The generalization that 
we consider is allowing for dipoles of arbitrary length. An 
"Z-dipole" is an object where the negative charge lies di- 
rectly / planes above its positive partner. While the orig- 
inal problem involved just the coupling of nearest neigh- 
bor planes, our generalized model involves all possible 
couplings. Associated with each Z-dipole is a fugacity 
yi/1-K (the 2ir is for convenience). The grand partition 
function for the system may be written as: 



{N k ,i} 



n 

k,l 



(N k ,iy. 



N, 



Z[{N k ,i}} 



(3.12) 



where N^i denotes the number of Z-dipoles in layer fc; N.i 
is the number of Z-dipoles in the system; and is the 
number of dipoles (of any length) that have their posi- 
tive ends in layer fc. The sum is over all particle number 
configurations {N/.j} that satisfy the charge neutrality 
constraint in each plane: = Y], Nk—u- The canon- 
ical partition function corresponding to a given dipolc 
distribution {Nf-j} is: 



Z{{N k ,i}} 



n 



'd 2 xt}d 2 x {2) 



(2) 



exp[-tf({iV M }) 



(3.13) 



H({Nk,i}) is the Hamiltonian (actually an entropy) cor- 
responding to the dipole distribution {Nk,i}- The coor- 
dinate Xu \ is the planar coordinate of the ith positive 

(2) 

charge of layer fc and x k \ is the planar coordinate of its 
negative partner which lives in layer fc + being 
the length of the dipole being described. The string con- 
straint is imposed by the delta function, where we use 
the normalization J R2 ^-S(f ) — 1. The product is over 
all positive charges in all layers. The integration is over 
the space il T . This is defined to be the set of all possible 
spatial configurations of the dipole distribution {Nk,i} 
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that respect the hard-core constraint: no two charges in 
a given plane may be closer than distance r. 

O ur p rocedure is an extension of the treatment in 
Refs. 17 20. The first part of an RG procedure normally 
involves integrating over the high momentum modes of 
the system. In our problem, these correspond to those 
configurations where in some plane, we have a pair of 
charges separated by a distance between r and r+dr. We 
assume a dilute system so only oppositely charged pairs 
are considered and also the distance between the mem- 
bers of a pair is taken to be much smaller than the dis- 
tance from the pair to another charge. The basic coarse- 
graining step in our RG transformation is illustrated in 

fig. El 



FIG. 5: The basic coarse-graining step in our RG transforma- 
tion. 

Suppose a particular state involves pairing the negative 
end of an Zi-dipole in layer k with the positive end of an 
Z 2 -dipole in layer k + 1%. Viewed at long length scales, we 
effectively have an Qi -H 2 )-dipole in layer k. We will find 
that integrating over all possible pairings gives a zeroeth 
order term (which just involves replacing f2 T with fl T+ dr) 
and a number of correction terms of order dr where two 
short dipoles were destroyed and replaced by a longer 
dipole. Since the procedure respects the charge neutral- 
ity constraint, these correction terms will combine with 
other terms in the grand partition sum. The second step 
involves rescaling lengths so that the high momentum 
cutoff, in the new variable, is the same as before. The 
aim is to see how the fugacities and couplings change as 
we run this procedure. 

Details of the calculation are given in Appendix B. 
Here we give the resulting flow equations: 



dyi 
dt 



(2 - n) yi 



dyi ,„ x 
dt 



l-l 

^ ' VmVl—rn 



m—1 



(3.14) 
(3.15) 
(3.16) 



where t = lnr. One notable feature is that the coupling 
does not change with the flow, in contrast with the 2d 



Coulomb gas where the coupling does vary (albeit at sec- 
ond order in the fugacity) . This indicates that strings are 
stiffer objects than charges. Another observation is that 
for the initial conditions of our physical problem, namely 
that yi(0) = y = 2Tre- E ^ kbT and yt(0) = for I > 1, 
the flow equations have an exact solution: 



yi = yoT 



1) 



l-l 



(3.17) 



Our RG is valid as long as the corrections to the fugacities 
are small, meaning that the derivatives dyi/dt should be 
bounded. If we look at the above result, Eq. 13.171 we 
see that when the term in brackets is greater than 1, yi 
diverges with I. Therefore, a critical length, which we 
interpret as an in-plane correlation length, is defined by 
when the term in brackets equals 1: 



RG 



1) = 1 



(3.18) 



Substituting earlier expressions and noting that for our 
system, k — 1/2, we find that: 



ln £|| ,RG 
£\\,RG 



1 



ln(e 



-E 3 /k B T 



K) 



9k B T V ' E s /k B T 
exp(32g / j lB JB/9k B T) (3.19) 



for the fields and temperatures of interest. This value 
is the same as that predicted in Ref. 11 using a free en- 
ergy argument. For r < iJG> yi decreases with I which 
means that states with long dipoles are less probable than 
states with short dipoles. If r > £||,-RGj Vi diverges with 
I which suggests that longer dipoles are favored, but, as 
mentioned above, the RG procedure is no longer valid 
in this regime. We note that when r = C||,HG> Vi is m ~ 
dependent of I so that, as in the mean field calculation 
discussed above, there is no discernible length scale in 
the z direction. 

If t < £ii rq, then we may consider an out-of-plane 
length scale, which we define nominally as the value of 
I = l T for which yi/yi = 1/e. 



l T = 1 



1 



/£' -1 

\ T 3 / 2 -l 



hi 



(3.20) 



We may interpret l T as the typical length of a string 
segment that is captured by a tube of diameter r (where 
a tube need not be straight). 



5. Comparison with simulation 

In Fig. we show the magnetization as a function of 
the magnetic field strength on a log-log scale. Our algo- 
rithm allows us to simulate spin ice in a [111] magnetic 
field with very high accuracy. 

The magnetization should scale with the average den- 
sity of defects, which in turn should scale like the inverse 
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FIG. 6: The crossover between exponents. 



square of the in-plane correlation length. As shown in 
this figure, the data at low fields are well fit by the expo- 
nent 8 /3 obtained in the mean field calculation discussed 
earlier. At somewhat higher fields, the data are well fit 
by the exponent 32/9, obtained by the RG calculation 
discussed earlier and also in Ref. [Tl] by looking at the en- 
tropic contribution to the free energy. At high fields, the 
exponent of 8L/3 (=16 for L=6 (sites), as was the case in 
the simulations) characterizes a regime where finite-size 
effects are important, as discussed below. 

The low field crossover makes qualitative sense in that 
at low fields, there will be many defects which screen 



one another which suggests that a mean field treatment 
may be reasonably accurate. At higher fields, the gas 
of defects is more dilute so an RG treatment would be 
required. 

The high field crossover is a finite-size effect since the 
position of a crossover between exponents is system size 
dependent and the corresponding exponent is also system 
size dependent, getting steeper with increasing system 
size. The finite-size behavior may be explained as fol- 
lows. At high magnetic fields, there are a small number 
of string defects in the system. The magnetization and 
the energy of one string defect in a system of size L are 
—4LgfiBJ/3 and ALg^BJB /3 respectively. The energy 
cost grows linearly with system size and, as mentioned 
above, the defects are favored solely due to their entropic 
contribution to the free energy. At sufficiently high mag- 
netic fields, a given system will be too small to provide 
the entropy to balance the energy cost of a string. This 
will occur when the magnetization per spin reaches the 
magnetization of a system with one string defect: 



= [1/3-2(4L/3)/(16L 3 )] 5 mb*/ 
= [l/3 — 1/(6L 2 )] gusJ- 



(3.21) 



In this case, the statistical weight of a sin- 
gle string defect will be a Boltzmann factor 
exp(—8LgiiBJB/3kBT) and the magnetization will 
equal [1/3 — C exp(— 8Lgfis JB/iksT)] g^sJ, where 
C is some constant. The crossover between different 
regimes occurs when the magnetization reaches l|3.21[l . 
We have good agreement with the 8L/3 behavior for a 
variety of system sizes, including L = 6 which is shown 
in figure H3 



IV. THE HIGH FIELD REGIME 

On the plateau, the magnetization of the triangular 
sublattice is saturated and we may consider each kagome 
plane separately. Thus, the 3-dimensional model may be 
mapped onto a 2-dimensional one. Whereas the spins 
in the triangular sublattice are fixed, the physics in the 
kagome planes remains non-trivial. Each triangle on the 
kagome plane contains two up pseudospins (cr = 1) and 
one down pseudospin (cr = —1). This Ising model on the 
kagome lattice may be mapped onto the dimer model on 
the hexagonal lattice ^LILH in which a down pseudospin 
corresponds to a dimer on the hexagonal lattice. The 
model retains an extensive ground state entropy, S/Ub = 
0.080765. 

If we flip a down (pseudo)spin it violates the ice rule. 
This corresponds to breaking a dimer into two monomers. 
As with string defects, these monomers may be separated 
and move freely on the lattice. The energy cost for cre- 
ating two monomers is 2£ = AJ c g — 2g^sJB /3. This 
energy vanishes at a critical field B c = 6J e s/(g^B J)- At 
higher fields the monomers proliferate leading to com- 
plete saturation and an ordered state with zero entropy. 
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The physics near the transition may be described by the 
following Hamiltonian which acts on the kagome lattice: 



H 
~T 



(4.1) 



where the sum is over all nearest neighbors; Si are classi- 
cal Ising spins taking values +1 and —1; his the strength 

of a fictitious magnetic field; and K++ = 0, = 

A'_+ = K = [ gf i B JB/6 - Jeff] /T, and K— = oo. The 
coupling constants imply that each triangle of the kagome 
lattice contains at most one down pseudospin and that 
down spins cost energy (positive or negative dependent 
on the magnetic field strength). 
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FIG. 7: The magnetization (top) and the entropy (bottom) 
around the transition between the plateaux. The simple 
Bethe approximation is compared to the Monte Carlo results. 
The exact result for the entropy at zero monomer density and 
Pauling's estimate for the entropy at zero magnetic field are 
shown for reference. The series expansion contains the results 
from Ref. ^| on the monomer-dimer model. 

We may calculate the magnetization and entropy using 
the simple Bethe approximation. Details are given in 
Refs. Il4ll2ll but we may quote the results: 

1 1 



2 1 + x 2 

3xz In z 
~2 + 6xz 



2z i 



iln- 

4 x 2 (3z 



x y 



(4.2) 
(4.3) 



where x = 2z/(l + \/l + 8z 2 ) and z = exp(-2K), 

In Fig. [3 we compare these expressions with a Monte 
Carlo simulation. The simulation is of a kagome lattice 
with 16x16 up-triangles (768 total spins). The standard 
single spin-flip Metropolis algorithm was used, which 
may explain the inaccuracy in the simulated entropy at 
low fields, where a more clever scheme may be needed to 
sample the degenerate manifold. The entropy was com- 
puted, for a given field, by integrating from high temper- 
atures (where = (3/4) ln2 per atom) to low tem- 
peratures. 

We find that the simple Bethe approximation is accu- 
rate for moderate and high monomer densities (higher 
fields) but does not work so well at low monomer den- 
sity (lower fields). As the Bethe approximation does not 
account for long cycles on the lattice, the approximation 
should indeed break down when the correlation length is 
large (monomer density is small). We note that the cor- 
relation length is infinite at zero monomer density since 
the dimer model on the hexagonal lattice is critical. 

In a higher-order series expansion, one may account for 
some corrections to the Bethe approximation^ As seen 
in the figure, the corrections are almost indiscernible for 
the magnetization. For the entropy, the corrections give 
better agreement at the low monomer density and are 
negligible at high monomer densities. 

There is a giant peak in the entropy at the transition 
point, S/ks = l/41n(16/5) « 0.291, which exceeds even 
the zero field entropy. The peak is due to the cross- 
ing of an extensive number of energy levels which have 
macroscopic entropies. For B = B c , the energies of states 
corresponding to different numbers of monomer defects 
are equal since the monomer and dimer weights are, by 
definition, equal at the critical field. There are an exten- 
sive number of states corresponding to a given number 
of monomers (below saturation) . The highly degenerate 
ground state manifold explains the large spike in the en- 
tropy. 



V. CROSSING POINTS 

The theory described in the previous section implies 
that the curves of magnetisation versus field, plotted for 
different temperatures, will display a crossing point. This 
arises simply because the partition function depends on 
magnetic field and temperature effectively only through 
the combination (B — B c )/T. Thus, when plotted as a 
function of B — B Cl the curves coincide only at the point 
B = B c . At this point, the Bethe approximation gives a 
value for the magnetisation of m = OAg^isJ, see Eg. 14.21 

In addition, we expect a crossing point at low fields, 
due the interplay of string and monomer defects. Indeed, 
where the plateau is well-formed, the string density is 
n s ~ exp(— 32g/is JB/9ki,T) and the monomer density is 
n m ~ exp(— 8f m /7fcsT), where £ — g^BJ{B c — B)/3 is 
the energy of creating one monomer. The crossing point 
occurs when n s = rib- With logarithmic accuracy, we can 



9 



write 

Z2gii B JB _ 8gfi B J(B - B c ) 

9k b T ~ 21k B T ' [ ' ' 

Thus the crossing point lies at B* = 3£? c /31. 

VI. RELATION TO EXPERIMENT AND 
OTHER THEORIES 

Our model gives a description of the high field transi- 
tion that is qualitatively consistent with experiment for 
a range of temperatures^. In particular, a peak in the 
entropy has been observed close to the high-field termi- 
nation of the plateau (Fig. 9 in Ref. |j). As this feature 
was taken to be an experimental artefact, it was not an- 
alyzed in detail in that work. However, it appears that 
its height is rather smaller than the one we find here, 
although the number of data points is not enough to de- 
termine the center of the peak or its height. 

However, recent experiments^ 2 on the spin ice com- 
pound Dy2Ti207 have indicated that at low tempera- 
tures, the high field transition becomes first order. In 
Ref. |22, the onset of first order behavior was found to 
occur for temperatures lower than a critical temperature 
of T c ~ 0.36K (~ 0.327J e ff, D y/k B )- Figure □ shows that 
our predicted curves remain continuous even at temper- 
atures below this observed T c . 

A likely reason for the discrepancy is the long range na- 
ture of the dipolar interaction, which we approximated 
as a nearest neighbor Ising model. The simplest way 
to account for these interactions is to model the ignored 
interaction terms as giving rise to a magnetic field pro- 
portional to the magnetization. By assuming the magne- 
tization M, as a function of the effective field B + aM, 
has the same functional form as given in figure[7| we may 
self-consistently determine M for a given B. Using a as 
a free parameter, we find that this simple model predicts 
the onset of first order behavior, at the experimentally 
observed critical field B c , only for temperatures in the 
millikelvin range. To obtain a higher numerical T c re- 
quires a larger a, which causes a lower numerical B c . To 
get the numerical T c to match experiment requires an 
a so large that our numerical B c is "negative" (in the 
sense of artificially extending the M = 1/3 line of figure 
0for the purpose of a spline fit). It seems that a more 
careful treatment of the dipolar interaction is required in 
order to explain the recent experimental results. Also, 
we have not considered the impact of the slowdown of 
the dynamics which is observed at low temperature^ 

As for the crossing points mentioned above, the high- 
field one does indeed appear to be present in the exper- 
imental data-* 2 ^ in the appropriate temperature range. 
The experimental value of the magnetization at the cross- 
ing point is about m = 0.38g/is J, reasonably close to the 
theoretical value m = 0Agfi B J. By contrast, a crossing 
point at small fields is harder to make out, and an ap- 
proximate estimate of its location gives B* = 0.35S C , in 



disagreement with the theoretical B* = 3£? c /31. 



VII. ENTROPY SPIKE AND 
MAGNETOCALORICS 

Fig. H shows a stark contrast between the behavior of 
magnetization and entropy as the field strength is in- 
creased. Whereas the magnetization increases monoton- 
ically going from one plateau to the other, the entropy 
displays a strong (but smooth) non-monotonicity. 

One question which naturally arises is whether such 
an entropy peak exists more generally between two mag- 
netization plateaus - what is the crucial ingredient for 
the existence of the spike? The sectors with different 
magnetizations are degenerate because not only do the 
monomer defects not cost any energy at the degeneracy 
point, but they also do not interact. Such a situation 
has in fact been observed already in a much more fa- 
miliar frustrated model, namely the triangular Ising an- 
tiferromagnet in a longitudinal field. Here, there is a 
(non-degenerate) low-field plateau with magnetization of 
1/3, in addition to the usual saturated high-field plateau. 
These two are separated by a degeneracy point where 
'up-up-up' and 'up-up-down' triangles are degenerated 
The statistical mechanics of that point is described by 
the hard-hexagon modelf 2 ^ the entropy of which is exten- 
sive. A similar phenomenon - a magnetization plateau 
bounded by two entropy spikes - also appears in the case 
of an effectively Id helimagnetiS& 

In classical Ising models, such a degeneracy seems not 
so surprising as the allowed energies are naturally dis- 
crete. However, a similar situation can arise even in 
bona-fide Heisenberg models. This follows from the re- 
sult by Richterei ai.jSl who demonstrated that near sat- 
uration, on a range of frustrated lattices (including the 
kagome), localized spin- 1/2 excitations exist. As one 
sweeps the magnetic field from saturation downwards, 
one would therefore also expect an entropy spike in those 
models. A numerical study testing this assertion is in 
progress^ 

Cooling by adiabatic (de)magnetization 

At low temperatures, near the degeneracy point, the 
partition function depends on magnetic field and tem- 
perature effectively only through the combination [B — 
B c )/T. One may thus argue that the spike may be used 
to effect cooling by adiabatic demagnetization 2 ^ in ex- 
actly the same way one may use paramagnets - analogous 
constraints limit the application in either case. 

There are two features which may be worth pointing 
out at this point. Both follow from the fact that - unlike 
in the case of a paramagnet - B c 7^ 0. Firstly, maximal 
cooling occurs at a finite field, namely around B c . This 
phenomenon may therefore be useful to effect cooling for 
a magnet in a field, with the restriction that B c , for a 
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given spin ice compound, is not tunable. Secondly, if B 
approaches B c from below, one can in fact obtain "cool- 
ing by adiabatic magnetization" , as entropy and magne- 
tization grow together in this regime. 
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VIII. CONCLUSIONS 

In this paper, we have analyzed in detail the magneti- 
zation curve of nearest-neighbor spin ice in a [111] mag- 
netic field. The basic ingredient which makes this sys- 
tem particularly interesting is that a uniform field can 
be used to couple to the Ising pseudospins 
gered fieldi22*2i This amounts to the possibility of apply- 
ing fields which would have appeared to be rather unnat- 
ural in the formulation of a simple Ising model (without 
the detour via spin ice) on the pyrochlore lattice. 

As a result, one observes an attractively rich behav- 
ior. Perhaps the most salient is the dimensional re- 
duction from pyrochlore to kagome under the appli- 
cation of an external field. The restoration of three- 
dimensionality upon weakening the field goes along with 
the one-dimensional string defects. We hope that the ex- 
tension developed here of Kosterlitz's RG treatment to 
such extended defects might be of more general use. 

A particularly attractive feature of the monomer- 
dimer model we have obtained here lies in the fact that 
the relative monomer and dimer fugacities in the low- 
temperature (T <C J c ff) regime are given by simple Boltz- 
man weights of Zeeman energies. They are thus straight- 
forwardly tunable by changing the strength of the ap- 
plied field. In particular, anisotropic fugacities can be 
obtained by tilting the field, and they therefore do not 
require an actual manipulation (such as an application 
of anisotropic stress) of the two-dimensional layer. 

As discussed previously in Ref^i the price for our 
ability to analyze the model in such detail has been 
the omission of the long-range nature of the dipolar in- 
teraction. A truncation of the interaction at only the 
nearest-neighbor distance would seem a rather drastic 
step; an expectation of quantitative agreement between 
experiment and the nearest-neighbor model will in gen- 
eral likely be misplaced. However, as we argue in a dif- 
ferent context, it turns out that, in an intermediate tem- 
perature regime, this is not entirely unreasonable. 13 This 
observation might lie at the basis of the fact that the mea- 
sured dipolar ice entropy agrees so well with Pauling's 
estimate. Our 'prediction' of the entropy peak between 
the intermediate and saturated plateaux bears witness 
to the promise of our approach to unearth at least some 
qualitative features of interest. 
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APPENDIX A: THE CLUSTER ALGORITHM 

We use a loop algorithm to simulate spin ice at low 
fields. The algorithm probes only the spin ice ground 
state manifold and therefore can work only at low tem- 
peratures and low magnetic fields. All attempted loop 
flips are accepted in our algorithm. 

The algorithm works as follows. To construct a loop, 
we, first, pick at random a tetrahedron of fixed orienta- 
tion (and mark it as a first tetrahedron in a loop), then 
we pick with probability 1/2 a spin direction (in or out 
of a tetrahedron) and pick a first spin in a loop using 
the following rules. If both spins with the chosen direc- 
tion are on the kagome sublattice then we pick the spin 
with a probability 1/2, which is independent of the spin 
orientation. If one spin is on the triangular sublattice 
and another is on the kagome sublattice then we pick 
the spin with probability that depends on the spin orien- 
tation. Namely, if the spin on the triangular sublattice 
is out of the tetrahedron (along the magnetic field), we 
pick the spin on the kagome or triangular sublattices with 
respective probabilities 



P2 = j±- (A2) 
1 + 9 

where g will be fixed by the detailed balance condition, 
see below, and p\ + p% = 1. If the spin on the triangular 
sublattice points into the tetrahedron, we pick the spin 
on the kagome or triangular sublattices with probabilities 
P2 and pi respectively. Then we flip the chosen spin thus 
introducing two defects in the tetrahedra that share the 
spin. 

After choosing the first spin, we move to the neigh- 
boring tetrahedron with a defect. The next tetrahedron 
has two spins with the opposite orientation. We flip one 
of these two spins adding it to the loop using the same 
prescription as we used to pick the first spin. Thus we 
move the defect to another tetrahedron. Then we repeat 
this procedure iteratively moving one of the two defects 
through the lattice until we encounter the other defect 
in the first tetrahedron - the two defects will annihilate 
and the loop will be closed. Since we add spins to the 
loop with alternating signs - two spins with opposite ori- 
entation from each tetrahedron we traverse, the ice rule 
is not violated. 

The algorithm is ergodic since any pair of different con- 
figurations differ by spins on closed loops only. They can 
always be connected by flipping these loops. 



11 



Let us sketch the proof of the detailed balance condi- 
tion. Suppose that we have flipped some loop. In order 
to prove detailed balance, the first site in a loop that re- 
turns us to the original configuration must be the first 
site in the original loop and the reversed loop must be 
constructed in the reverse direction. We can prove the 
detailed balance condition locally, i.e. for all short se- 
quences of the loop, see Fig. |S] and Fig. |5J It is easy to 
check that most of these sequences are trivial, i.e. they 
have equal energies before and after spin flip and equal 
probabilities to go from one to another configuration. An 





P(A — > B') = pi/2 and the reverse probability of going 
from B' to A is P(A <- B') = p 2 /2. We have from fOjl 



g = — = exp(— 8/i/3). 
Pi 

Therefore if we choose p\ and p 2 as 

1 



and 



Pi 



P2 



1 



-8/i/3 ' 



,-8h/3 



1 + e -Sh/3 

then the detailed balance condition is fulfilled. 

APPENDIX B: RG CALCULATION 



(A4) 



(A5) 



(A6) 



FIG. 8: Configurations A and B. Tetrahedra are shown on 
top of each other. Small arrows indicate a short sequence of a 
loop. Up and down spins are denoted by black and grey dots. 

example of such a simple sequence is shown in Fig. |H1 
The probability of going from configuration A to config- 
uration B is equal to the probability of going from B to 
A (equal to 1/2). In order to prove the detailed balance 
condition, we only need to consider the energies of sin- 
gle spins that are the second spins in the sequences (the 
energies of the first spins in the sequences are taken into 
account in the previous step). These spins have the same 
energies. Thus the detailed balance condition is satisfied 
trivially. An example of a nontrivial sequence is shown 
in Fig. [5] The energies of configurations A and B' are 
different there. We have to prove the detailed balance 



We introduce the abbreviation: 





FIG. 9: Configurations A and B' . Tetrahedra are shown on 
top of each other. Small arrows indicate a short sequence of a 
loop. Up and down spins are denoted by black and grey dots. 

condition 

P(A -> B) /P{A <— B') = P{B')/P{A) (A3) 

The right hand side in l|A3(l is just a ratio of Boltz- 
mann weights and is equal to exp(8/i/3), where h — 
gUsJB /HbT ', since the energy of configuration A (the 
energy of the second and third spins in the sequence) is 
4hksT/3, and the energy of configuration B' (the en- 
ergy of the second and third spins in the sequence) is 
—AhksT '/3. According to our algorithm, the probabil- 
ity of going from configuration A to configuration B' is 



(1) d 2 x 



(2) 




(1) (2) 



(Bl) 



in terms of which the canonical partition function for 
a given dipolar distribution {Nk,i} may be written: 
Z({N kt i},T) = f n &l T exp(-H). Our RG calculation 
has two steps. The first step is integrating over short 
length scales, i.e. those states where at least one pair of 
charges is separated by a distance between r and r + dr. 
The second step is to rescale variables to restore the short 
distance cutoff. When we carry out the first step, the re- 
sult is a zcroth order term and a correction of order dr: 



Z({N kil },r) 



d£l T exp(— H) 



klmij 



(B2) 

where Ikimij is the contribution of the configuration that 
has the negative end of the ith m-dipole of layer k paired 
with the positive end of the jth (I — m)-dipole of layer 
k + m. The sum over k is over all planes; the sum over I 
is over all dipole lengths up to the number of planes; and 
the sum over m is from 1 to I — 1. The form of this term 
is given by: 



l-klmij 



n' 



d(x {2) ,t) 



d^ ( xf-xf 



(B3) 

The region of integration of the positive charge is 
an annulus of radius t and thickness dr centered on the 
negative charge x\ . This region is denoted by d{x\ ,t). 
The position of this negative charge (and hence the pair) 
is integrated over the entire area A. Strictly speaking, 
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(2) 

x\ ' would have to avoid the hard cores of all of the other 
charges but this introduces an error of order (dr) 2 . Q T+dT 
is the space of configurations of the rest of the charges 
in which the charges are separated from each other by 



(2) 



.r : ') refers to the 



a distance of at least r + dr. H(x 

(2) 

piece of the Hamiltonian which involves charges x\ and 
and the rest of the Hamiltonian is denoted by H' . 
The x^ integration amounts to making the substitu- 



tion Xj — x) 



d 2 x^ = rdrdO; and integrating 



over angles. If we denote the latter two of integrals of 
equation IB3I by /, then: 



(B7) 



where the space ^ t j_'^™ is analogous to f2 T +dTj except that 
there is one less m-dipole in layer k; one less (I — m)- 
dipole in layer k + m; and one more ^-dipole in layer k. 
What we are actually interested in is the grand partition 
function feouation I3.12|l . Because our RG procedure is 
consistent with the charge neutrality constraint, the var- 
ious {Ikimij} may be combined with different terms in 
the grand partition function. When we substitute into 
eg uat ion 13 . 1 21 and arrange terms, we find that: 



I = — 

T 



dr f dPx? ] ^ 



(1) (2) 



(B4) 



We assume that our gas of defects is sufficiently dilute 
that the following distances are much greater than the 
pair separation r: (1) the distance of a particle in plane 
k+m from our pair, (2) the distance of a particle in plane 

k from the positive charge x+ , and (3) the distance of a 

(2) 

particle in plane k + I from the negative charge Xj . In 
this dilute limit, we may make the approximation: 



/r (1) t (2 \ rx {2) 



?(2) 



-(1) -(2) 
X\ ' — X - ' 



(Wo) 



We also have that H(x[ 2 \ Xj) is small in this limit, 
which allows us to expand the exponential and to leading 
order, the integral may be done exactly 17 . The result is: 



/ = —61— 

r V r 



,(2) 



(2vr 



Ee a e b In — 
T 



(B6) 



In the penultimate line, the sum refers to a sum over all 
charges, positive and negative, residing in the plane k + 
to. This sum term may be neglected in the large A limit, 
which is why, in contrast to the Kosterlitz calculation 20 , 
the coupling strength does not vary during our RG flow 
(see equation 13. 14|) . The delta function implies that the 
TO-dipole and (I — m)-dipole have been combined into a 
larger Z-dipole. Returning to our correction term: 



n(tr 



dft T exp(-H) 



k,l,m k,\ 



2TT 



dr y m yi 
t (2tt) 



N h ,-1 



(B8) 



The prime on the second product means that y l 
has been taken outside the product. If the fugacities are 
small, then we may write this in a more convenient way: 



(yi + t- Emii ymVi- m ) N - 1 



e [n 

{N kA } k,l 
/ dfi T exp(— H) 



{2k) n *.i(N ki )\ 



(B9) 



Finally, we rescale lengths, x —* x(l + dr/r) , and find 
(dropping primes): 



Z = 



[n 



ft) ' 



dO T exp(-i/) (BIO) 



where 



i-i 



i , dr v—v v / di~ . . dr . ,„ 

Vi = (yi + — V y m yi- m ){i + 2—1 - k—) Bli 



'klmij 



2tt— 
r 



fc.i,: 



exp(-iJ) 



The flow equations l|3.14|l follow from this. 
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